load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "~/work/mylib/callvar.ncl" 

begin

  tp = "all" 
  tp = "tropic" 
  tp = "polar" 
  tp = "midlat" 

  load "load_tp" 

  fna = "K2009_"+tp+"_Ni.dat" 
  fnb = "K2009_"+tp+"_T.dat" 

  ni = asciiread(fna,-1,"float") 
  st = asciiread(fnb,-1,"float")

  nt = dimsizes(ni) 

  ni!0 = "time" 
  ni&time = ispan(1,nt,1) 

  st!0 = "time" 
  st&time = ispan(1,nt,1) 

  
  ni_210 = where(st.gt.205. .and. st.lt.215, ni, -999) 
  ni_210@_FillValue = -999

  ni_220 = where(st.gt.215. .and. st.lt.225, ni, -999) 
  ni_220@_FillValue = -999

  ni_230 = where(st.gt.225. .and. st.lt.235, ni, -999) 
  ni_230@_FillValue = -999

  ni_240 = where(st.gt.235. .and. st.lt.245, ni, -999) 
  ni_240@_FillValue = -999

  n210 = vnmiss_1D(ni_210) 
  n220 = vnmiss_1D(ni_220) 
  n230 = vnmiss_1D(ni_230) 
  n240 = vnmiss_1D(ni_240) 

;;...........................................................................
;; 210K 
;;...........................................................................
  qsort(n210)

  dm = dimsizes(n210)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1

  pv = n210 
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 205K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pa = pp 

;;...........................................................................
;; 220K 
;;...........................................................................
  qsort(n220)

  dm = dimsizes(n220)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n220
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 215K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pb = pp 

;;...........................................................................
;; 230K 
;;...........................................................................
  qsort(n230)

  dm = dimsizes(n230)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n230
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 225K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pc = pp 

;;...........................................................................
;; 240K 
;;...........................................................................
  qsort(n240)

  dm = dimsizes(n240)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n240
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 235K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pd = pp 

;;...........................................................................
;; output data 
;;...........................................................................
  fno = "K2009_"+tp+".nc" 

  system(" rm "+fno) 

  flo = addfile(fno,"c") 

  flo->ni=ni  
  flo->st=st
  

;;...........................................................................
;; output data 
;;...........................................................................
  yval = new((/4,6/),"float",-999.)

  do is = 0,5
     yval(0,is) = pa(is)
     yval(1,is) = pb(is)
     yval(2,is) = pc(is)
     yval(3,is) = pd(is)
  end do 
 
  x = (/210., 220., 230., 240./)		


;;...........................................................................
;; wks
;;...........................................................................
  wks = gsn_open_wks("pdf","box_"+tp)

;;...........................................................................
;; res
;;...........................................................................
  res            = True

  res@tiMainString = ""

  res@tmXBLabels = (/"205K-215K","215K-225K","225K-235K","235K-245K"/) ; labels for each box
  res@tmXBLabelFontHeightF = 0.016

;;...........................................................................
;; resources for polylines that draws the boxes
;;...........................................................................
  llres                   = True			
  llres@gsLineThicknessF  = 2.5                 ; line thickness 

;;...........................................................................
;; resources that control color and width of boxes
;;...........................................................................
  opti          = True			
  opti@boxWidth = 1.6
  opti@boxColors = (/"blue","red","green","black"/) ; Color of box(es)
  ;;res@gsnMaximize = True                        ; Maximize box plot in frame.

  plot = boxplot(wks,x,yval,opti,res,llres)	; All 3 options used...



;;...........................................................................
;; add some polymarkers
;;...........................................................................
  mres               = True                     ; marker mods desired
  mres@gsMarkerIndex = 1                        ; polymarker style
  mres@gsMarkerSizeF = 50.                      ; polymarker size
  mres@gsMarkerColor = "black"                    ; polymarker color
  xi  = yval(0,5)                               ; polymarker locations
  xi2 = yval(1,5)
  xi3 = yval(2,5)
  xi4 = yval(3,5)

  dum1 = gsn_add_polymarker(wks,plot,210.,xi,mres) 
  dum2 = gsn_add_polymarker(wks,plot,220.,xi2,mres) 
  dum3 = gsn_add_polymarker(wks,plot,230.,xi3,mres) 
  dum4 = gsn_add_polymarker(wks,plot,240.,xi4,mres) 
  
  mres@gsMarkerColor = "navy blue"              ; change color
  mres@gsMarkerIndex = 1                        ; change style
  mres@gsMarkerSizeF = 40.                      ; change size

;;  dum4 = gsn_add_polymarker(wks,plot,-3.,xi-4.,mres) 
;;  dum5 = gsn_add_polymarker(wks,plot,-1.,xi2-2.,mres)
;;  dum6 = gsn_add_polymarker(wks,plot,1.,xi3-1.,mres) 
  draw(plot)
  frame(wks)
  
end




 
